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I.  INTRODUCTION 


The  cylinder  expansion  test,^  in  which  the  explosively-driven,  outward 
radial  expansion  of  a  standard  metal  cylinder  is  observed  by  a  streak  camera, 
has  become  one  of  the  classic  experimental  tools  in  research  concerned  with 
detonation  dynamics.  In  fact  a  large  portion  of  our  practical  and  reliable 
information  about  explosives  has  been  obtained  directly  or  indirectly  from 
this  experiment;  Gurney  velocities,  Taylor  angles  and  detonation  velocities 
can  all  be  deduced  from  cylinder  expansion  test  data.  In  practice,  however, 
the  primary  use  of  such  data  has  been  in  determining  equations  of  state  (EOS) 
for  the  detonation  product  gases,  expressed  as  a  curve  in  the  pressure-density 
plane  describing  the  Chapman-Jouguet  isentrope.  This  is  normally  accomplished 
by  an  iterative  procedure  in  which  a  full-scale  continuum  mechanics  computer 
code,  such  as  HEMP,  is  used  to  model  the  observed  experiment.  The  equation 
of  state  in  the  codes  is  written  in  an  assumed  analytic  expression 
involving  a  number  of  undetermined  parameters  which  can  be  adjusted  repeatedly 
with  successive  computer  calculations  until  "good"  agreement  with  the  exper¬ 
imental  measurements  is  obtained.  The  resulting  parameter  values  are  then 

considered  as  correct  for  the  oarticular  exnlosive  a™1  can  reasonably  be  used 
in  subsequent  calculations.  Viewed  in  this  light, the  computer  codes  can  be 
seen  as  very  elaborate  fitting  and  extrapolation  procedures. 

Among  the  better  known  expressions  for  the  equations  of  state  of  gases 
are  the  polytropic  gas  law  (constant  y) ,  LJD  (Lennard-Jones-Devonshire) ,  BKW 
(Becker-Kistiakowsky-Wilson) ,  and  JWL  (Jones-Wilkins-Lee) .  Most  of  the 
presently  available  codes  will  allow  the  user  to  select  the  form  which  he 
considers  most  appropriate  and  also  contain  a  compilation  of  appropriate 
parameter  values  for  commonly  used  explosives.  (See  References  2  and  3.) 

W.  Kury  et  al .,  "Metal  Acceleration  by  Chemical  Explosives,"  Proceedings 
of  the  Forth  Symposium  on  Detonation,  Office  of  Naval  Research  Report 
ACR-126,  pp3-12  (1965). 

2 

Dobratz,  Brigitta  M.,  "Properties  of  Chemical  Explosives  and  Explosive  Stim¬ 
ulants,"  UCRL  51319 ,  Dec  1972,  Lawrence  Livermore  Laboratory ,  Livermore  CA. 

3 

Lee,  E.L.,  Homig,  H.C.,  Kury,  J.W.,  "Adiabatic  Expansion  of  High  Explosive 
Detonation  Products,  "  UCRL  50422,  May  1968,  Lawrence  Livermore  Laboratory, 
Livermore,  CA. 
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In  the  present  study  we  propose  an  alternative  approach  for  deducing 
equation  of  state  data  directly  from  the  cylinder  expansion  test  results. 

This  procedure  is  basically  a  reversed  version  of  the  well-known  analysis  of 

4 

cylindrical  bombs  due  to  G.  I.  Taylor.  It  does  not  require  the  use  of  a 
large  scale  continuum  mechanics  code,  nor  does  it  assume  an  analytic  expres¬ 
sion  for  the  equation  of  state  with  adjustable  parameters,  but  computes 
pressure  directly  in  tabular  form  as  a  function  of  density.  To  implement 
this  procedure  a  computer  program  called  •'CEDAR*'  (Cylinder  Expansion  Data 
Reduction)  was  written  in  the  BASIC  programming  language  for  the  Hewlett 
Packard  Model  9845B  minicomputer;  execution  time  for  this  program  is  typi¬ 
cally  less  than  one  minute.  A  tabulation  of  CEDAR  is  contained  in  Appendix 
A  and  an  example  of  an  output  is  shown  in  Appendix  B. 

II.  EXPERIMENTAL  SET  UP 

To  provide  a  general  framework  for  discussion,  let  us  first  describe  the 
experimental  set  up  and  procedure  for  the  cylinder  expansion  test  as  it  is 
conducted  in  the  Warhead  Mechanics  Branch  of  the  Ballistics  Research  Labor¬ 
atory.  A  schematic  is  shown  in  Figure  1. 

A  long,  hollow, open-ended  cylinder  is  packed  with  explosive  to  a  desired 
loading  density  and  a  detonator  and  booster  charge  placed  at  one  end.  The 
most  common  choice  for  this  test  is  a  12  inch  long  copper  cylinder,  with  a 
one  inch  interior  diameter  and  0.1  inch  wall  thickness.  Interior  diameters 
of  two  and  four  inches  are  also  considered  standard  and  results  from  one 
geometry  can  be  made  comparable  with  others  through  appropriate  similarity 
variables,  as  in  Reference  1.  The  program  listed  in  Appendix  A  can  be  used 
for  any  choice  of  inner  and  outer  cylinder  diameters.  In  the  tests  performed 
at  Lawrence  Livermore  National  Laboratories  great  care  was  taken  to  fabricate 
the  cylinders  with  exact  dimensions.  At  BRL, much  more  economical  stock 

*G.I.  Taylor ,  " Analysis  of  the  Explosion  of  a  Long  Cylindrical  Bomb  Detonated 
at  one  End,  "  The  Scientific  Papers  of  Sir  G.  I.  Taylor,  Volume  III,  p277 ; 
edited  by  G.K.  Batchelor,  Cambridge  University  Press,  1963. 


copper  pipe  was  purchased  from  local  sources  and  machined  to  proper  dimensions, 
practically  identical  results  were  obtained. 

A  Beckman  and  Whitley  Model  168  streak  camera  is  placed  behind  a  blast 
wall  and  situated  to  observe,  using  a  mirror,  the  transverse  expansion  of  a 
chosen  cross  section  of  the  cylinder;  this  section  is  located  so  that  a  fully 
developed  steady-state  detonation  front,  with  no  end  effects,  prevails  during 
the  period  of  observation.  For  the  standard  12  in.  cylinder, the  observed 
station  was  positioned  9.S  in.  from  the  initiated  end. 

As  the  detonation  wave  passes  through  the  observed  section, the  cylinder 
walls  are  driven  outwards  and  obscure  the  backlighting  provided  by  an  argon 
light  bomb.  The  external  diameter  at  that  section  is  continuously  observed 
by  the  camera  through  a  slit  aperture  and  its  image  swept  at  uniform  velocity 
across  the  recording  film  by  mean;  of  a  rotating  drum  mirror.  A  typical  film 
negative  resulting  from  this  procedure  is  shown  in  Figure  2.  The  vertical 
scale  is  proportional  to  the  cylinder  wall  displacement.  (Both  top  and  bottom 
walls  are  observed  and  have  virtually  identical  behavior.)  The  proper  scale 
factor  for  converting  film  measurements  to  actual  displacements  is  determined 
from  the  millimeter  scale,  appearing  at  the  left,  which  was  photographed  before 
initiation.  Measurements  along  the  horizontal  axis  correspond  to  elapsed  time; 
the  appropriate  scale  factor  can  be  determined  from  knowledge  of  the  camera 
dimensions  and  mirror  rotation  rate. 

When  a  film  record,  such  as  that  in  Figure  2,  becomes  available  from  an 
experiment fit  is  placed  in  a  film  reader  and  x  and  y  coordinates  are  read  from 
both  top  and  bottom  displacement  curves.  These  are  compared  for  discrepancies, 
averaged  and  multiplied  by  scale  factors  to  yield  a  curve  d(t)  which  expresses 
the  exterior  cylinder  wall  displacement  as  a  function  of  time. 

As  an  additional  feature,  the  detonation  velocity  of  the  explosive  can 
be  determined  during  the  experiment  by  placing  several  electrical  contact  pins 
or  break  wires  at  measured  distances  along  the  cylinder.  At  our  facility  four 
contact  pins  were  positioned  at  2  in.  intervals  beginning  at  2.5  in.  from  the 
initiating  end.  The  elapsed  time  between  the  making  or  breaking  of  contact 


at  different  stations  is  recorded  as  voltage  on  an  analog  device  during  the 
experiment  and  quantified  afterwards  using  an  oscilloscope.  Sane  care 
must  be  taken  to  ensure  that  the  explosive  is  packed  at  uniform  density 
along  the  entire  length  of  the  pipe.  In  our  case  this  was  accomplished  by 
forming  the  explosive  separately  in  2  in.  sections  with  the  proper  loading 
density.  These  were  then  placed  firmly  together  in  the  cylinder,  but  without 
additional  compaction. 


III.  ANALYSIS 

Determination  of  an  equation  of  state  for  the  detonation  product  gases 

requires  an  understanding  of  the  explosion  gas  dynamics,  its  interaction  with 

the  cylinder  walls  and  the  material  strength  of  the  cylinder  at  high  strain 

4 

rates.  The  discussion  by  G.  I.  Taylor  provides  a  useful  and  workable  model 
of  this  process  and  will  serve  as  a  basis  for  the  present  analysis.  In 
addition,  we  have  developed  a  methodology  for  extracting  the  cylinder  wall 
motion  in  more  detail  than  is  directly  available  from  the  experimental  record. 
Together  these  enable  us  to  attain  our  basic  objective  —  determination  of 
pressure  as  a  function  of  density  (or  specific  volume)  along  the  Chapman- 
Jouguet  isentrope.  Our  approach  involves  three  subtasks  which  will  be  dis¬ 
cussed  separately  in  the  following  sections;  kinematics  of  the  cylinder  wall 
motion,  parametric  determination  of  the  EOS  and  considerations  relating  to 
numerical  analysis. 

A.  Kinematics  of  the  Cylinder  Wall  Motion 

The  curve  d  ■  d(t)  obtained  from  the  experimental  film  readings  is  a 
record  of  the  apparent  cylinder  wall  displacement  and  not  a  direct  measurement 
of  the  motion  of  an  individual  wall  particle.  This  measured  displacement  can  be 
considered  as  Eulerian  whereas  the  path  of  an  actual  wall  particle  is  more 
properly  described  by  a  Lagrangian  coordinate  system  in  which  both  the  axial 
and  radial  displacements  are  expressed  as  functions  of  time.  In  order  to 
compire  wall  particle  trajectories  directly  with  code  predictions  and  to 
compute  accelerations  more  accurately , it  is  necessary  to  deduce  this  Lagrangian 
information  from  knowledge  of  d(t)  only.  In  the  following  discussion  we  shall 
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show  how  this  can  be  done. 

To  begin,  we  define  our  terms  more  precisely.  Let  z  and  r  denote  axial 
and  radial  coordinates  in  an  Eulerian  system  chosen  so  that  z  =  0  at  the 
observation  station  and  r  =  0  along  the  central  axis  of  the  cylinder. 

We  suppose  that  the  inner  and  exterior  radii  of  the  test  cylinder  at  station 
z  *  0  are  given  by  R  (t)  and  R  (t)  which  are  initially  Cat  time  t  £  0)  equal 
to  Rq  and  Rj  respectively.  The  experimental  data  measures  the  displacement 
of  the  outer  surface  from  rest  ;  thus> 


Ra  (t)  -  Rj  +  d(t) 


To  determine  R  (t)  we  assume  incompressibility  of  the  cylinder  wall  material 
£1 

during  acceleration.  This  implies 

u(Re  "  r2J  =  Ao  5  *(R?  * 
e  a  o  l  o 


-  [rJ  (t)  -  A0/»]  1/2  . 


The  acceleration  of  the  wall  should  be  calculated  not  at  the  inner  or 
outer  surface  but  at  the  central  surface  defined  by 


5  fRe"’  -  V2"] 


Let  us  next  assume  that  the  expansion  is  steady  state,  that  is  it 
propagates  with  a  fixed  waveform  along  the  cylinder  at  a  constant  velocity  D, 
the  explosive  detonation  velocity.  Thus  the  motion  at  any  axial  station  z  is 
related  to  that  at  z  >*  0  through  the  similarity  variable  (t  -  z/D)  ;  in 
particular,  the  apparent  central  surface  displacement  from  its  initial  position 
at  station  z  observed  at  time  t  is 


h(z,t)  = 


)(t-z/D)-Rm  CO) 


t  £  z/D 
t  >  z/D . 


1 


r 


A.)  as  viewed  by 


B.)  as  viewed  by 
Streak  Camera 
at  station  2»  0 


r 
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icle  Motion,  Viewed  in  Eulerian 


B.) 


Figure  5.  Comparison  of  Two  Different  Particle  Trajectories 

a)  particle  observed  at  z  =  0  at  time  t  b)  particle 
originating  at  z  =  0,  observed  at  time  t'  =  i,(t)/D 


Considered  as  a  function  of  z  only, this  describes  the  central  wall  dis¬ 
placement  as  it  would  appear  to  a  framing  camera  activated  at  time  t;  as  a 
function  of  t,it  gives  the  displacement  that  would  be  viewed  by  a  streak 
camera  located  at  station  z.  These  two  different  frames  of  reference  are 
shown  schematically  in  Figure  3. 


The  wall  particle  initially  located  at  station  z  =  0  is  driven  outwards 
along  a  path  which  might  be  described  parametrically  in  the  form 


z  =  f(t) 

r  =  g(t)  +  Rm(0),  (4) 

as  indicated  in  Figure  4.  The  steady  state  assumption  establishes  a 
relationship  between  these  functions  when  d(t)  is  known  since 


g(t)  =  h(f(t),t)  =  Rm(t  -  f (t)/D)  -  Rm(0)  (5) 


However,  this  is  insufficient  to  determine  f  and  g  uniquely  unless  an 

4 

additional  assumption  is  made.  Accordingly ,we  suppose t  as  in  Taylor’s 
analysis,  that  pressure  is  exerted  normally  and  consequently  that  no 
longitudinal  deformation  of  the  cylinder  walls  occurs.  Mathematically  this 
means  that  the  distance  between  wall  particles  measured  as  arc-length  along 
the  central  surface  remains  constant  during  the  expansion.  We  therefore 
introduce  the  function 


*(t) 


[1  ♦  (— ) 
3? 


1/2 

]  dz 


which,  referring  to  Figure  5. a,  is  seen  as  the  arc-length  observed  at  time 
t  between  the  detonation  front,  at  z  =  Dt,  and  the  observation  station 
z  =  0.  After  a  change  of  variables  this  can  be  more  conveniently  written  as 


The  implication  of  our  last  assumption  is  that  the  wall  particle  observed 
passing  through  z  =  0,  r  =  R^ft)  at  time  t  must  have  originated  at  station 
z  =  Dt  -  Jl(t)  on  the  cylinder  wall.  If  the  cylinder  is  viewed  moments  later 
at  time. 


t*  =  «,( t)/D  >  t  for  t  >  0  t 

we  see  that  the  expansion  has  progressed  to  the  situation  indicated  in  Figure  5.b. 
From  this  figure  we  can  now  deduce  the  relationship 

f(f)  =  f(£Ct)/D)  =  £(t)  -  Dt  (7) 

and  also,  by  comparison  with  Figure  5. a, 

g(t')  =  h(0,t).  (8) 

Moreover,  the  physical  significance  of  t'becanes  apparent  -  it  is  the  time  elapsed 
since  the  particle  observed  at  station  z  =  0,  at  time  t  first  commenced  its  motion. 
This  might  be  considered  as  time  in  a  Lagrangian  system. 

For  any  choice  of  t  we  may  compute  P(t)  and  t'  from  (6)  and  knowledge  of 
d(t)  thereby  obtaining  from  (7)  and  (8)  the  position  (f(t'),g(t')  ♦  Rm(0))  which 
is  reached  by  the  cylinder  wall  particle  initially  located  at  z  «  0,  after  time 
t'  following  passage  of  the  detonation  wave.  By  doing  this  for  sufficiently  many 
choices  of  t  we  can  reconstruct  the  entire  particle  path  in  the  z,r  (Eulerian) 
coordinate  system.  Henceforth  we  consider  motion  in  the  Lagrangian  system,  using 
t'  as  the  independent  variable  of  time;  t  can  be  treated  as  a  parameter. 

B.  Gas  Dynamics  and  Parametric  Determination  of  the  Equation  of  State 

Having  obtained  a  complete  history  of  the  cylinder  wall  motion  the  next 
step  is  to  reconstruct  the  Chapman-Jouguet  (C-J)  isentrope  for  the  explosive 
detonation  products.  This  is  a  curve  P  =  P(V)  ,  relating  pressure  to  specific 
volume,  which  describes  the  adiabatic  expansion  from  the  C-J  state  (Pcj>vcj) 
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outwards  as  far  as  the  experimental  data  permits.  Normally  pressure  decreases 
from  several  hundred  to  just  a  few  kilobars.  Briefly,  our  procedure  for 
accomplishing  this  is  as  follows:  1)  differentiation  of  the  particle  dis¬ 
placements  to  obtain  velocity  and  acceleration  histories;  2)  calculation 
of  pressure  as  a  function  of  acceleration  and  cylinder  wall  strength;  3) 
determination  of  density  as  a  function  of  expansion  radius  and  gas  velocity. 
The  apparent  radial  expansion  Re>  measured  at  discrete  choices  of  time,  t,  is 

used  parametrically  in  this  process,  while  R  ,  R  ,  t',  f(t'),  g(t'),  P  and  V 

m  a 

are  determined  subsequently  for  each  choice  of  Rg(t),  thereby  providing  one 
point  (P,V)  of  the  isentrope. 

The  first  step,  differentiation  of  the  displacement  data,  requires 
special  numerical  treatment  and  will  be  discussed  separately  in  the 
following  section. 

The  pressure  of  the  detonation  product  gases  at  any  moment  follows 
directly  from 


Force  _  Ma 
Area  2nl^ 


(9) 


where  M  is  the  mass  of  the  cylinder  per  unit  length  and  a  is  the  instantane¬ 
ous  acceleration,  measured  at  the  center  of  mass  of  a  wall  particle.  This 
is  essentially  the  same  as  the  expression  derived  by  Taylor  from  a  different 
point  of  view.  The  material  strength  of  the  cylinder  walls  has  been  ignored 
here  but  could  be  taken  into  account  by  introducing  additional  terms;  for 
example, 

,o  -  R  ,  R  -  R 

P  =  Ma  +  Ll - 1)  tk  -  kin  -£ - S-  )  (10) 

2wRa  2,*  m  5  Rl-R0 
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where  and  kt  denote  the  engineering  and  true  linit  strength  of  the  wall 
Material.  These  additional  teras  should  really  be  considered  as  empirical 
adjustments  since  little  is  known  about  the  strength  of  materials  at  strain 
rates  of  10^  or  10^  per  second.  For  standard  copper  cylinders  assumed  values 
of  k^  or  ks  ,less  than  10  kilobars  have  only  minor  effects  on  the  pressure/ 
density  calculations  which  we  carried  out; even  at  50  kilobars ,the  influence 
of  the  additional  terms  in  (1)  is  only  seen  at  the  lower  end  of  the  P(V) 
curve  (P  <  5  kilobars) .  Such  values  are  far  beyond  the  measured  strength  of 

copper  at  lesser  strain  rates  so  we  may  assume  that  material  strength  is  not 
a  critical  consideration  at  this  point. 

After  pressure  has  been  determined  as  a  function  of  the  radial  displacement, 

4 

Rg,  the  specific  volume  follows  from  equations  (3)  and  (11)  of  Taylor  ;  in  our 
notation  these  become 


V  =  -  = 

P 


<y 


where  p  is  the  explosive  loading  density,  u  is  the  gas  velocity  given  by 


u  =  D  + 


2TrpQR^D 


© 


pf(f) 

2  .  i  2 

+  [  dg(t’)l  1  ; 

(13) 

L  dt' 

L  dt'  J 

and  w  is  the  cylinder  wall  velocity  measured  at  the  central  surface 


w2(t) 


the  functions  f(t')  and  g(t')  are  determined  as  described  above. 

In  carrying  out  this  pressure/density  calculation  it  is  helpful  but 
not  necessary  to  know  the  detonation  pressure  in  advance  either  from 
a  separate  experiment  or  analytical  predictions,  since  this  provides  a 
starting  point  for  plotting  in  the  P-V  plane.  Such  a  priori  knowledge  is 
usually  assumed  in  the  parameter  fitting  method  of  EOS  formulation,  as 
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discussed  earlier.  If  PCJ  is  known, then  the  specific  volume  at  detonation  can 
be  determined  from  (11)  and  (12)  with  w  =  0  and  R  =  R ;  thus, 

a  O 


(DPr 


A  useful  quantity  which  can  also  be  computed  after  the  other  variables 
have  been  determined  is  the  generalized  ratio  of  specific  heats 


r  =  - 


(9lnP\ 

alnV/c 


This  can  be  approximated  numerically  by 


r.  = 

i 


KV  pi_i  ] 
^vi.i/vi3 


where  i  is  an  index  of  consecutive  points  on  the  C-J  isentrope.  However, 
this  quantity  is  a  bit  delicate  to  compute  since  it  implies  third  order 
differentiation  of  the  displacement  data  and  can  be  expected  to  produce 
eccentric  results.  Nevertheless,  with  the  data  modifications  outlined  in 
the  next  section,  we  have  been  able  to  achieve  a  fairly  decent  resolution 
of  this  function.  In  contrast,  it  must  be  said,  an  analytic  expression  of  the 
EOS,  with  fitting  parameters,  makes  T  easy  to  evaluate  and  smoothly  behaved. 

A  good  deal  of  commentary  has  been  directed  at  interpreting  the  morphology 
of  r  in  physical  terms,  but  this  mostly  seems  to  miss  the  point,  since 
the  appearance  and  behavior  of  r  is  as  much  a  consequence  of  the  assumed 
EOS  form  as  it  is  of  any  underlying  physical  properties.  From  mathematical 
analysis  we  know  that  two  functions  can  closely  agree  and  still  have  very 
dissimilar  derivatives.  Thus, it  is  possible  to  add  terms  to  the  EOS  formula 
which  cause  little  effect  on  subsequent  displacement  calculations  but  which 
strongly  affect  the  appearance  of  r  . 


C.  Numerical  Differentiation  and  Data  Adjustment 

As  anyone  who  has  attempted  it  knows,  differentiation  of  data  obtained 
from  experimental  measurements  is  frought  with  difficulties,  double  differen¬ 
tiation  invites  disaster  and  third  order  derivatives  are  usually  chaotic. 

This  would  appear  to  be  true  in  the  present  case  even  though  the  basic  input 
data  appears  as  smooth  as  might  reasonably  be  hoped  for. 

If  discrete  measurements  t.  and  d(t^)  are  taken  from  the  experimental  film 
record  and  indexed  by  i,  then  let  Rm  ^  and  tj  denote  the  resulting  values  of  Rm 
and  t'  determined  as  explained  in  Section  A.  The  usual  numerical  approximations 
for  the  radial  components  of  the  velocity  and  acceleration,  derived  from  (4)  and 
(8) ,  are 

R  -  R  » 

v.  =  -li-iti - (17) 

1  t'  -  t’. 

l+l  l 


a.  =  2  ^ - Vi~1-  ■  (18) 

1  +*  +* 

z  i+1  "  z  i-1 

with  analogous  formulas  for  axial  motion. 

These  expressions  are  exact  if  the  displacement  is  quadratic  in  time,  which 
generally  is  not  the  case.  When  these  formulas  were  used  in  conjunction  with 
(10)  and  (11)  to  compute  the  C-J  isentrope  for  Comp  B,  grade  A,  using  the  test 
data  published  in  Reference  I  (See  Table  l)fwe  obtained  the  results  indicated 
in  Figure  6;  the  JWL  equation  of  state  curve  is  shown  for  comparison. 

Since  this  was  of  little  use  we  had  either  to  abandon  or  modify  the 

basic  procedure.  Fortunately,  through  numerical  experimentation,  we  were 
able  to  obtain  greatly  improved  results  with  only  minor  corrections  to  the 
input  data  and  with  a  plausible  modification  of  the  differentiation  formulas. 


The  first  point  to  be  realized  is  that  the  published  data,  shown  in  columns 
1  and  2  of  Table  1,  lacks  precision.  The  displacement  time  measurements  listed 
in  column  2  are  given  in  units  of  .01  microseconds.  Although  this  may  represent 
the  greatest  accuracy  achievable  with  the  available  equipment  it  implies  round¬ 
off  errors  of  +  .005  microseconds,  not  to  mention  the  inaccuracies  in  measurement. 
This  may  not  be  significant  in  itself  but  does  exert  considerable 
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influence  on  the  higher  order  derivatives.  It  seems  reasonable  therefore,  and 
not  a  violation  of  scientific  procedure,  to  introduce  small  adjustments 
(£  .005  microseconds)  if  this  improves  the  coherence  of  the  subsequent  cal¬ 
culations;  After  all,  physical  insight  tells  us  that  the  acceleration  of  the 
walls  should  be  a  reasonably  smooth  function,  monotonically  decreasing  with 
time  and  displacement.  Calculations  which  show  a  widely  oscillating  acceleration 
must  therefore  be  spurious. 


TABLE  1.  ORIGINAL  AND  ADJUSTED  DATA  FROM  LLNL  COPPER  CYLINDER  TEST 

REFERENCE  1. 


R-RO 

(mm) 


TIME 

(micro  sec) 


MODIFIED  TIME 
(micro  sec) 


2 

3 

4 

5 

6 


2.17 

3.00 

3.77 

4.51 

5.22 


2.17 

2.997 

3.76845 

4.5045 

5.2160 


7 

8 
9 

10 

11 


5.91 

5.90957 

6.59 

6.5895 

7.26 

7.25875 

7.92 

7.91941 

8.57 

8.57265 

12 

13 

14 

15 

16 


9.22 

9.21929 

9.86 

9.95998 

10.50 

10.49527 

11.13 

11.12567 

11.75 

11.7516 

17 

18 

19 

20 
21 


12.37 

12.3736 

12.99 

12.9921 

13.60 

13.60735 

14.22 

14.21958 

14.83 

14.8290 

22 

23 


15.43  15.43575 

16.04  16.04 
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TABLE  2.  CALCULATED  VALUES  FOR  ACCELERATION 
AND  T  DURING  CYLINDER  EXPANSION 


Expm ' t 

Displ cm't 

(1) 

Original  Data; 
formulas  (1716(18) 

Modified 

formulas 

(2) 

Data ; 
(17)6(18) 

(3) 

Modified  Data; 
formulas  (21)5(22) 

R-Ro 

Acc 

r 

Acc 

r 

Acc 

r 

0.0 

.000 

2.49 

.000 

2.49 

.000 

2.49 

1.0 

.773 

3.07 

.773 

3.07 

.720 

2.41 

2.0 

.252 

3.94 

.252 

3.94 

.261 

3.60 

3.0 

.141 

3.51 

.146 

3.35 

.145 

3.55 

4.0 

.109 

2.19 

.100 

2.89 

.099 

2.93 

5.0 

.060 

4.48 

.073 

2.71 

.072 

2.78 

6.0 

.070 

-.67 

.056 

2.66 

.054 

2.71 

7.0 

.049 

3.32 

.043 

2.66 

.042 

2.70 

8.0 

.023 

7.10 

.034 

2.62 

.033 

2.65 

9.0 

.025 

-.36 

.027 

2.62 

.027 

2.63 

10.0 

.027 

-.31 

.023 

2.55 

.022 

2.55 

11.0 

.029 

-.31 

.020 

1.79 

.020 

1.73 

12.0 

.005 

20.22 

.018 

1.48 

.018 

1.47 

13.0 

.032 

-27.81 

.017 

1.38 

.017 

1.38 

14.0 

.004 

24.68 

.016 

1.36 

.016 

1.36 

15.0 

.034 

-35.93 

.015 

1.48 

.015 

1.49 

16.0 

.036 

-.30 

.014 

1.40 

.014 

1.40 

17.0 

.003 

31.42 

.012 

2.19 

.012 

2.22 

18.0 

.003 

2.18 

.011 

2.06 

.011 

2.03 

19.0 

.039 

-63.55 

.011 

1.33 

.010 

1.28 

20.0 

.04-3 

-1.44 

.010 

1.38 

.010 

1.38 

21.0 

.039 

2.23 

.010 

1.04 

.010 

1.03 

22.0 

.041 

-.46 

.009 

.99 

.009 

1.03 

23.0 

— 

— 

— 

— 

.009 

.93 

These  calculations  have  been  done  three  different  ways  as  noted:  Either 
the  original  or  modified  displacement  data  from  Table  1  was  used  t  and 
either  formulas  (17)  and  (18)  or  formulas  (21)  and  (22)  were  used  for  the 
acceleration  calculation. 
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After  some  trial  and  error  adjustments  of  the  original  data, we  obtained  the 
modified  data  indicated  in  column  3  of  Table  1  as  a  more  precise  record  of  the 
actual  cylinder  wall  displacement.  The  adjusted  and  unadjusted  values  of 
acceleration  versus  time  are  given  in  Table  2.  It  is  clear  here  that  the 
modified  values  exhibit  a  much  more  reasonable  behavior  than  the  original 
ones  and  are  probably  a  more  accurate  description  of  the  actual  physical 
event.  The  resulting  P-V  curve  is  shown  in  Figure  7.  We  wish  to  stress  that 
the  modified  data  in  column  3  does  not  contradict  the  original  data  but  merely 
enhances  it  -  makes  it  more  precise.  Except  in  three  cases, the  adjustment  is 
less  than  the  round-off  error  of  .005  microseconds  and  the  maximum  adjustment 
is  only  .007  microseconds.  Essentially, we  have  used  the  data  to  correct 
itself  in  a  manner  which  leads  to  monotonically  decreasing  second  derivatives. 

Since  a  fair  amount  of  discrepancy  still  remained  between  our  calculated 
values  and  the  JWL  curve,  as  shown  in  Figure  7,  we  considered  whether  additional 
changes  in  the  procedure  might  yield  improved  results.  One  possibility  was 
to  change  the  velocity  and  acceleration  formulas  (17)  and  (18);  as  mentioned, 
these  are  exact  for  displacements  which  are  quadratic  in  time,  implying  a 
constant  acceleration  over  the  time  interval.  In  actuality , the  acceleration  is 
monotonically  decreasing  with  respect  to  time  and  displacement,  so  we  might 
suppose  that  over  a  time  interval  t  <  t  £  tj,  the  displacement  is  better 
described  in  the  form 

Rm(t)  =  do  +  Vo(t  -  V  +  C(t  -  Vb  (19) 

where  dQ  and  vq  are  the  displacement  and  velocity  at  t  *  t  .  Positive,  strictly 
decreasing  acceleration  requires  that 

1  <  b  <  2  and  c  >  0 
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Figure  7.  C-J  Isentrope  Calculated  Using  Adjusted  Data  and  Algorithm 
Based  on  Equations  (17)  and  (18). 
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and  in  fact  b  =  1.5  seems  to  produce  the  most  satisfactory  results  in  practice. 

T.  ,  .  are  known  and  if 

If  d  and  vo 
o 

W  *  di 

is  specified  then 


c  =  (dj  -  dQ  -  vQ  At)/ (At) 


(20) 


where  At  »  tj  -  tQ.  The  acceleration  at  mid-interval,  t  =  (tQ  +  tj)/2,  is 


a  =  22-b(b-l) (dj-dQ  -  vQAt)/At2  (21) 

and  the  velocity  at  t  =  tj  can  then  be  evaluated  as 

Vj  =  vQ  +  aAt.  (22) 

At  the  start  of  the  first  time  interval f the  displacement  and  velocity  are 
known  to  be  zero;  the  average  acceleration  a  and  the  velocity  Vj  can  thus 
be  determined  in  terms  of  dj.  Formulas  (21)  and  (22)  can  then  be  applied 
recursively  to  determine  the  average  acceleration  over  each  subsequent  time 
interval . 

When  this  revised  procedure  is  applied  to  the  modified  displacement 
data  from  Table  1 ,we  obtain  the  acceleration  history  indicated  in  Table  2 
and  the  EOS  curve  shown  in  Figure  8.  This  appears  to  be  a  significant 


PRESSURE 


Figure 


.01  1 _ I _ I - 1 - 1 - 1 - 1 - 1 - 1 - 1 

.5  1  2  4  8 

RELATIVE  VOLUME  CV/V0) 

8.  C-J  Xsentrope  Calculated  Using  Adjusted  Data  and  Algorithm 
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improvement  over  the  more  straightforward  approach  using  equation  (17)  and  (18) 
(assuming  the  JWL  equation  of  state  for  Comp  B,  grade  A,  is  valid) .  We 
believe  that  it  is  noteworthy  that  no  further  modification  of  the  displacement 
data  was  necessitated  by  this  fundamental  change  in  algorithm  and  take  this 
as  additional  evidence  of  the  suitability  of  the  modified  data. 


IV.  DISCUSSION  AND  CONCLUSIONS 


The  procedure  described  in  this  report  provides  a  means  for  extracting 
equation  of  state  data  directly  from  cylinder  expansion  test  measurements 
and  should  provide  a  useful  tool  to  researchers  in  explosion  dynamics.  Its 
principal  advantage  is  in  eliminating  the  dependence  on  large  scale  continuum 
mechanics  code  calculations  and  contrived  mathematical  expressions  for  the 
equation  of  state.  It  must  be  said  howeve^that  this  procedure  still  is  a 
long  way  from  being  confidently  established  and  additional  work  is  required 
to  accomplish  this. 

Our  objective  in  the  present  report  has  been  to  derive  and  elucidate  the 
mathematical  formulations  and  methodology.  Thus  we  have  only  considered  the 
single  case  of  Comp  B,  grade  A,  in  a  standard  copper  cylinder.  To  complete 
the  task  a  thorough  comparison  of  this  procedure  against  known  results  needs 
to  be  carried  out.  In  particular: 

1)  The  method  should  be  applied  to  the  cylinder  test  data 
for  a  variety  of  explosives  and  results  compared  with 
conventional  EOS  data. 
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2)  The  EOS  subroutines  of  continuum  mechanics  codes  should 
be  modified  to  accept  the  tabular  data  produced  by  the 
present  method,  in  place  of  the  usual  analytic  expression; 
calculations  could  then  be  performed  and  compared  with  the 
previous  results  and  with  experiments. 

3)  Sensitivity  of  the  procedure  to  variation  in  experimental 
parameters  (cylinder  material  and  dimensions,  explosive 
loading  density)  should  be  investigated. 

4)  A  method  for  modifying  the  basic  input  data  should  be 
systematized,  if  possible,  so  that  experimenters  do  not 
have  to  waste  time  on  trial-and-error  data  adjustments.  In 
this  regard  it  may  be  possible  to  predetermine  the  properties 
of  r  (see  equation  (15));  that  is,  to  carry  out  data  adjust¬ 
ments  subject  to  restrictions  on  the  range,  fluctuation  or 
asymptotic  behavior  of  T. 

From  the  mathematical  point  of  view  we  believe  that  there  is  an 
interesting  and  practical  question,  worthy  of  a  separate  investigation, 
which  underlies  the  development  of  the  present  algorithm  and  relates  to 
the  general  problem  of  data  differentiation.  We  might  state  this  as 


follows: 


Let  the  values  of  a  function  be  given  at  discrete  choices 
of  its  argument,  let  E  >  0  be  a  given  error  bound  and 
N  >  0  be  an  integer.  How  ,  and  when,  can  the  given  data 
be  approximated  by  a  fitting  function  that  differs  from 
the  prescribed  values  by  at  most  E  and  which  possesses 
a  non-negative  derivative  of  order  N  in  the  region  of 
interest? 

The  data  adjustment  procedure  of  Section  III.C  amount  to  an  empirical 
solution  of  this  problem  but  the  question  is  certainly  of  broader 
interest . 
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TABULATION  OF  COMPUTER  PROGRAM  "CEDAR" 


TABULATION  OF  COMPUTER  PROGRAM  "CEDAR" 


1  !  RE-STORE "CEDAR :F" 

2  PRINT  “**********************************************************************" 

3  PRINT  M****************** ***********  CEDAR  »★★★★****★* **************** 

4  PRINT  *****************  CYLINDER  EXPANSION  DATA  REDUCTION  PROGRAM  ********** 

5  PRINT  "***************************************************************♦*****♦" 

6  PRINT 

7  PRINTER  IS  16 

8  !  PRINTER  IS  0 

10  REM  Unit  of  length:  millimeters 

20  REM  Unit  of  time:  microsecond 

30  REM  Unit  of  mass:  gram 

40  REM  Unit  of  velocity:  mm/mlcrosec  ( -Km/sec) 

50  REM  Unit  of  acceleration:  mm/mlcrosec*? 

60  REM  Unit  of  density:  gm/mm®3 

70  REM  Unit  of  pressure:  terapascal  (-  10  megabars) 

80  REM  Input  pressure  In  glgapascals  (=  10  kllobars) 

90  REM  Output  pressure  displayed  In  glgapascals 

100  DIM  D1 sp( 100) ,Ta( 100) ,Ca( 100 ) ,Rm( 100) ,Tp( 100) ,Gamna( 100) ,Ps( 100) 

110  DIM  T(100),Ra(100),Z(100},V(100),P(100).Rho(100) 

120  DIM  Vr(100),Vz(100),Vgas(100),Vo1(100) 

130  DIM  Re( 100) ,Ve( 100) ,Ae( 100) ,Acc(100) 

200  D-7.98  !  DETONATION  VELOCITY 

210  Dc-. 001717  !  EXPLOSIVE  LOADING  DENSITY 

220  Dm-. 008940  !  CYLINDER  WALL  DENSITY 

240  Km-1  !  WALL  STRENGTH  PARAMETER  (EQUATION  10) 

250  Ks*l  !  WALL  STRENGTH  PARAMETER  (EQUATION  10) 

260  REM  ****  ENTER  Km.Ks  IN  GIGAPASCALS  **** 

270  Psc-.OOl 

280  Km*Km*Psc 

290  Ks»Ks*Psc 

300  REM  ****  Km.Ks  NOW  IN  COMPATIBLE  UNITS  (terapascal s)  ********************* 
310  RO-12.7 

320  Rl-RO+2.606 

330  N-26 

340  PO-29.5  !  EXPLOSIVE  DETONATION  PRESSURE 

350  REM  ****  ENTER  PO  IN  GIGAPASCALS  **** 

360  PO-PO*Psc 

370  REM  ****  PO  NOW  IN  COMPATIBLE  UNITS  (terapascal s)  ************************ 

380  AO-PI*( Rl®2-RO®2 )  !  CROSSECTIONAL  AREA  OF  CYLINDER  WALL 

390  RmO«SQR( (R0*2+Rl®2)/2)  !  INITIAL  RADIAL  COORDINATE  OF  CENTRAL  SURFACE 

400  M«Dm*AO  !  CYLINDER  MASS  PER  UNIT  LENGTH 

410  C-Dc*PI*R0®2  !  EXPLOSIVE  MASS  PER  UNIT  LENGTH 


Note:  The  symbol  "*"  In  this  tabulation  denotes  exponentiation. 
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510  FOR  1*1  TO  N 

520  READ  01 sp(I),Ta(I) 

530  Re(I)=Rl+D1sp(I) 

540  NEXT  I 

550  DATA  0.00, -2. 00, 0.00, -1.00,0.00, -.140. 1.00, 1.212 

551  DATA  2., 2. 17 

552  DATA  3. ,3.00 

553  DATA  4., 3. 77 

554  DATA  5., 4. 51 

555  DATA  6., 5. 22 

556  DATA  7., 5. 91 

557  DATA  8., 6. 59 

558  DATA  9., 7. 26 

559  DATA  10., 7. 92 

560  DATA  11., 8. 57 

561  DATA  12., 9.22 

562  DATA  13., 9.86 

563  DATA  14., 10. 50 

564  DATA  15., 11. 13 

565  DATA  16., 11. 75 

566  DATA  17., 12. 37 

567  DATA  18., 12. 99 

568  DATA  19., 13. 60 

569  DATA  20., 14. 22 

570  DATA  21., 14. 83 

571  DATA  22., 15.43 

572  DATA  23., 16.04 

600  REM  **********  INPUT  DATA  SMOOTHING  FACTORS  ************** 

601  !  GOTO  700  !  SKIP  TO  700  TO  OMIT  DATA  SMOOTHING 

610  FOR  1*1  TO  N 

620  READ  Ca(I) 

630  Ta(I)*Ta(I)+Ca(I) 

640  NEXT  I 

650  DATA  0,0,0,0,-0.0040 

651  DATA  -.00300, -.00155, -.00550, -.00400, -.00043 

652  DATA  -.00050,-. 00125, -.00059, +.00265, -.00071 

653  DATA  -.00002, -.00473, -.00433,+.00160,+.00360 

654  DATA  +.00210,+.00735,-.00042,-.00105,+. 00560 

655  DATA  -.00035 

700  PRINT  "*********  INPUT  CYLINDER  EXPANSION  DATA  *********" 

710  PRINT  "  I  R-RO  T  Va  Aa" 

720  Ve(0)«0 

730  Ae(l)*0 

740  Rm(l)*RmO 

750  Ra(l)*RO 

800  REM  **********  CALCULATE  APPARENT  VELOCITY  AND  ACCELERATION  ************** 

810  FOR  1*1  TO  N 

820  Ve(I)*(Re(I+l)-Re(I))/(Ta(I+l)-Ta(I))  !  EXTERIOR  NALL  VELOCITY 

830  Ae(I)«2*(Ve(I)-Ve(I-l))/(Ta(I+l)-Ta(I-l))  !  EXTERIOR  NALL  ACCELERATION 

840  Rm(I)«SQR(Re(I)«2-A0/(2*PI))  !  CENTRAL  SURFACE  DISPLACEMNT 

850  Ra(I)*SQR(Re(I)C2-A0/PI)  i  INTERIOR  NALL  DISPLACEMENT 

860  PRINT  USING  861;I,D1sp(l).Ta(I),VeU),Ae(I) 

861  IMAGE  DDD,4(0DD.DDDD) 

870  NEXT  I 

880  PRINT 
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890  PRINT  "  I  *  INDEX  OF  DATA  POINT  READING  FROM  FILM  RECORD" 

900  PRINT  "R-RO  *  MEASUREO  EXTERIOR  WALL  DISLACEMENT  (MM)" 

910  PRINT  *  T  *  MEASUREO  EXPANSION  TIME  (MICROSECONDS)" 

911  PRINT  "  Va  *  APPARENT  EXTERIOR  WALL  VELOCITY  (MM/MICROSEC)" 

912  PRINT  "  Aa  =  APPARENT  EXTERIOR  WALL  ACCELERATION  0*/MICR0SEC®2)" 

913  PRINT 

1000  REM  ************  CALCULATE  LAGRANGIAN  TIME  SCALE  Tp( I )  ******************* 

1001  PRINT  "*************************  CALCULATED  DATA  *********************** 
1010  Tp(l)=0 

1020  Z(1)*0 

1030  Vr(l)=Vz(l)*V{l)*0 

1040  P(1)=P0 

1050  Vgas(l)*0 

1060  Vol (l)*l/Dc-P0/(D*Dc)®2 

1070  Rho(lM/Vol(l) 

1080  Gamma  ( 1 )  »Dc*D*>2/P0- 1 

1090  PRINT  "  I  Rm-RmO  Z  Tp  Vr  Vz  Ar  Ac  Rho  P 

Gamma ” 

1091  IMAGE  DDD,3(DDD.DDD),2X,3(DDD.DDD),DDD.DDD,D.DDDDD,DDDD.DD,DDD.DD 

1092  GOTO  1200  !  SKIP  TO  1200  TO  USE  IMPROVED  ACCELERATION  ALGORITHM 

1100  REM  *********  CALCULATION  BASED  ON  STANDARD  ACCELERATION  ALGORITHM  ******* 
1110  Tdel  *SQR(  ( Ta ( 2 ) -Ta(  1 )  )®2+(  ( Rm( 2 ) -Rm(  1 ) ) /D )®2 )  !  LAGRANGIAN  TIME  INCREMENT 

1112  Tp(2)*Tp(l )+Tdel  !  TIME  ELAPSED  IN  MOTION 

1114  Z(2)*D*(Tp(2)-Ta(2)+Ta(l))  !  LAGRANGIAN  AXIAL  DISPLMNT 

1120  FOR  1*2  TO  N-l 

1122  Tdel =SQR( (Ta(I+l)-Ta(I) )®2+( ( Rm( 1+1 ) -Rm( I ) ) /D )r2 )  !  LAGRANGIAN  TIME  INCR. 


!  EXPLOSIVE  PRODUCTS  PRESSURE 


1124  Tp(I+l)*Tp(I)+Tdel  !  TIME  ELAPSED  IN  MOTION 

1126  Z(I+l)*D*(Tp(I+l)-Ta(I+l)+Ta(l))  !  LAGRANGIAN  AXIAL  DISPLACEMENT 

1130  Vr( I )*{Rm( 1+1 ) -Rm( I ) )/Tdel  !  LAGRANGIAN  RADIAL  VELOCITY 

1132  Vz(I)*(Z(I+l)-Z(I))/Tdel  !  LAGRANGIAN  RADIAL  VELOCITY 

1134  V(I)=SQR(Vr(I)®2+Vz(I))  !  VELOCITY  ALONG  PARTICLE  PATH 

1140  Ar=2*(Vr(I)-Vr(I-l))/(Tp(I+l)-Tp(I-U)  !  LAGRANGIAN  RADIAL  ACCELERATION 

1142  Az*2*(Vz( I )-Vz( 1-1) )/(Tp( I+1)-Tp( 1-1))  !  LAGRANGIAN  AXIAL  ACCELERATION 

1144  Acc(I)*SQR(Ar®2+Az®2)  !  ACCELERATION  ALONG  PARTCLE  PATH 

1146  IF  Ra(I)<=Ra(I-l)  THEN  1180 
1150  P(I)=M*Acc(I)/(2*PI*Ra(I)) 

1152  P( I  )=P(I )*(1+(V( I )/D)®2)*( -1.5)  !  EXPLOSIVE  PRODUCTS  PRESSURE 

1154  Rde1*Re(I)-Ra(I) 

1156  P( I)*P(I)-Rde)*(Km-Ks*LOG(Rdel /<  Rl-RO) ) )/{ 2*PI*Ra( I ) ) 

1160  Vgas( I )*D+M*V( I )«2/(2*PI*D*Dc*R0«'2)-P( I)*(Ra( I )/R0)®2/(D*Dc) 

1162  !  Vol (I)=l+.5*(M/C)*(V(I)/D)®2-P(I)/(Dc*D®2)*(Ra( I)/R0)®2 

1164  !  Vol  ( I )=Vol ( I )/Dc*(Ra( I )/R0)®2 

1166  Vol { I )*Vgas( I )*Ra( I )r2/( Dc*D*R0®2 )  !  SPECIFIC  VOLUME 

1168  Rho( I )*l/Vol (I)  !  GAS  DENSITY 

1170  IF  P(I)>0  THEN  1174 
1172  P( I)*. 000001 

1174  Gamma( I  )*LOG(P( I )/P(  1-1 )) /LOG (Vol  (I-D/Vol  (I))  !  GENRL.  RAT. OF  SPEC. HEAT 

1176  GOTO  1190 
1180  P(I)*P(I-1) 

1182  Vaas(I)*Vgas(I>l) 

1184  Vol ( I )*Vol ( 1-1) 

1186  Rho(I)*l/Vol(I) 

1188  Gamma(I)*Ganma(I*-l) 

1190  PRINT  USING  1091;I,Rm(I).Rm0tZ(I),Tp(I),Vr(I),Vz(I),Ar,Az,Rho(I),P(I)/Psc, 
Gamma(I ) 

1192  NEXT  I 
1194  GOTO  1500 
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1200 

1205 

1210 

1220 

1230 

1240 

1242 

1250 

1252 

1260 

1270 

1272 

1280 

1282 

1290 

1300 

1310 

1312 

1314 

1320 

1322 

1324 

1330 

1340 

1342 


REM  *********  CALCULATION  BASED  ON  IMPROVED  ACCELERATION  ALGORITHM  ******* 
FOR  1-2  TO  N 

Tdel-SQR( (Ta{ I )-Ta( 1-1 ) )«'2+( (Rm( I )-Rw{ 1-1) )/D)®2)  !  LAGRANGIAN  TIME  INCR. 


Tp(I)*Tp(I-l)+Tdel  ! 

Z(I)«D*(Tp(I)-Ta(I)+Ta(l))  « 

Ar=Rm{ I)-Rm(I-l) -Vr( I -1 )*Tde1 
Ar*1.06*Ar/Tde1®2  ! 

Az»Z(I)-Z(I-l)-Vz(I-l)*Tdel 
Az=1.06*Az/Tdel®2  ! 

Acc( I )«SQR(Ar®2+Az®2)  ! 

Vr( I )-Vr( I-l)+Ar*Tde1  ! 

Vrm=Vr(I-l)+.5*Ar*Tdel  ! 

Vz( I )-Vz( I-l)+Az*Tdel  ! 

Vzm-Vz(I-l)+.5*Az*Tdel  ! 

V(  I)«SQR(Vrn£>2+Vzia®2)  ! 

IF  Ra(I)<-Ra( 1-1)  THEN  1400 


TIME  ELAPSED  SINCE  MOTION  COMMENCED 
LAGRANGIAN  RADIAL  ACCELERATION 

LAGRANGIAN  RADIAL  ACCELERATION 

LAGRANGIAN  AXIAL  ACCELERATION 
ACCELERATION  ALONG  PARTICLE  PATH 
RADIAL  VELOCITY,  INDEX  I 
RADIAL  VELOCITY  AT  MIDPOINT 
AXIAL  VELOCITY,  INDEX  I 
AXIAL  VELOCITY  AT  MIDPOINT 
AVERAGE  VELOCITY  ALONG  PATH 


Raa=.5*(Ra( I-1)+Ra( I))  ! 

Rea».5*(Re(I-l)+Re(I))  ! 

Rdel*Rea-Raa  ! 

P(I)=M*Acc(I)/(2*PI*Raa)  ! 

P(I)-P{I)*(1+(V(I)/D)®2)®(-1.5) 

P( I)=P( I ) -Rdel * ( Km- K  s*LOG ( Rdel /(Rl-R0)))/(2*PI*Raa) 

Vgas( I )-D+M*V( I )®2/( 2*PI*D*Dc*R0®2 ) -P( I )*{ Raa/R0)®2/{ D*Dc) 
!  Vol ( I )»1+.5*(M/C)*( V( I )/D)®2-P( I )/(Dc*D®2)*(Raa/R0)®2 
!  Vol ( I )=Vo1 ( I )/Dc*(Raa/R0)®2 


DISPLACEMENT 

DISPLACEMENT 


AVERAGE  INNER  HALL 
AVERAGE  OUTER  NALL 
WALL  THICKNESS 
DETONATION  PRODUCTS  PRESSURE 


1346  Vol (I)=Vgas(I)*Raa®2/(Dc*D*R0®2)  I  GAS  SPECIFIC  VOLUME 

1350  Rho(I)*l/Yol(I)  !  GAS  DENSITY 

1360  IF  P(I)>0  THEN  1380 
1370  P( I)-. 000001 

1380  Gamma( I)-LOG(P(I)/P( I-l))/LOG(Vol (I-D/Vol (I))  !  GENRL .RATIO  OF  SPEC. HE AT 
1390  GOTO  1450 
1400  P(I)«P(I-1) 

1410  Vga$(I)*Vgas(I-l) 

1420  Vol(I)»Vol(I-l) 

1430  Rho(I)-l/Vol(I) 

1440  Gamnta(I)-Gainna(I-l) 

1450  PRINT  USING  1091;I,Rm(I)-Rm0,Z(I),Tp(I),Vr(I),Yz{I)tAr,Az.Rho(I),P(I)/Psc, 
Gamma(I) 


1460 

1500 

1510 

1520 

1530 

1540 

1550 

1560 

1570 

1580 

1590 

1600 

1610 

1620 

3000 

3010 

3011 

3012 


NEXT  I 
PRINT 

PRINT  "Rm-RwO 
PRINT  "  Z 


PRINT 

N 

Tp 

PRINT 

N 

Vr 

PRINT 

W 

Vz 

PRINT 

N 

Ar 

PRINT 

M 

Az 

PRINT 

N 

Rho 

PRINT  "  P 
PRINT  "  Ganna 
PRINT 

PRINT  "END  OF 


»  RADIAL  DISPLACEMENT  OF  HALL  PARTICLE  CENTER  OF  MASS  (PM)" 

-  AXIAL  DISPLACEMENT  OF  NALL  PARTICLE  CENTER  OF  MASS  (MM)" 

-  ELAPSED  TIME  SINCE  START  OF  PARTICLE  MOTION  (MICROSEC)" 

-  LAGRANGIAN  RADIAL  VELOCITY  COMPONENT  (W/MICROSEC)" 

-  LAGRANGIAN  AXIAL  VELOCITY  COMPONENT  (PW/MICROSEC)" 

-  LAGRANGIAN  RADIAL  ACCELERATION  COMPONENT  (W/MICR0SEC®2)" 

-  LAGRANGIAN  AXIAL  ACCELERATION  COMPONENT  (m/MICR0SEC®2)" 

-  DETONATION  PRODUCTS  GAS  DENSITY  ( GRAMS />*®3)" 

-  DETONATION  PRODUCTS  GAS  PRESSURE  (GIGAPASCALS)” 

-  GENERALIZED  RATIO  OF  SPECIFIC  HEATS" 

COMPUTATION" 


REM  ***************  PLOT  PRESSURE  VS  SPECIFIC  VOLUME  ********************* 
GRAPHICS 

!  GOTO  3020  !  SKIP  TO  3020  TO  USE  9872A  PLOTTER 

PLOTTER  IS  13 ."GRAPHICS" 
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3013  LIMIT  0,180,0,131 

3014  GOTO  3030 

3020  PLOTTER  IS  7, 5. "9872 A" 

3021  LIMIT  0,200,25,250 
3030  LOCATE  25,95,15,80 
3040  CSIZE  3 

3050  SCALE  LOG( .5) ,L0G(32) ,LOG( .01) ,L0G(100) 

3060  CLIP  LOG( .5) ,LOG( 16) ,LOG( .01) ,LOG( 100) 

3070  AXES  L0G(SQR(2)),L0G(10000)/8,L0G(.5),L0G(.01) 

3071  OUTPUT  705;"SM*“ 

3072  LINE  TYPE  1 

3073  FOR  1-2  TO  N-l 

3074  PLOT  L0G(Dc*Vo1 (I)),LOG(P(I)/Psc)  !  PLOT  DATA  POINT  ALONG  C-J  ISENTROPE 

3075  NEXT  I 

3080  OUTPUT  705; "SH" 

3084  PENUP 

3085  LINE  TYPE  1 

3090  REM  PLOT  C-J  ISENTROPE  FROM  JWL  EQUATION  OF  STATE 

3091  A-5. 24229 

3092  B-. 076783 

3093  Rl-4.2 

3094  R2-1.1 

3095  Omega-. 34 

3096  FOR  1-1  TO  N-l 

3097  Ps( I )-A*EXP( -Rl*Vol { I )*Dc)+B*EXP( -R2*Vo1 ( I )*Dc)+.0108*( Yol ( I )«Uc)®( -(1+Ome 
ga)) 

3098  PLOT  L0G(VoI(I)*0c),L0G(Ps(I)*100) 

3099  NEXT  I 

3100  FOR  1-0  TO  5 

3110  MOVE  L0G(2®( 1-1 ) )+LOG( .9) ,LOG( .006) 

3120  LABEL  2®(I-1) 

3130  NEXT  I 

3140  FOR  1-0  TO  4 

3150  MOVE  L0G(.3),(I-2)*L0G(10) 

3160  LABEL  10*(I-2) 

3170  NEXT  I 

3180  MOVE  L0G(1),L0G(.003) 

3190  LABEL  "RELATIVE  VOLUME  (V/VO)" 

3200  LOIR  PI/2 

3210  MOVE  LOG( .25), LOG! .1) 

3220  LABEL  "PRESSURE  (GlgaPascal s)" 

3221  LOIR  0 

3230  MOVE  L0G(3),L0G(10) 

3231  LABEL  "J.H.L." 

3232  MOVE  L0G(3) ,L0G(5) 

3233  LABEL  "EQUATION" 

3234  MOVE  L0G(3),L0G(2.5) 

3235  LABEL  "OF  STATE" 

3250  MOVE  L0G(.7),L0G(.4) 

3251  LABEL  "REOUCEO” 

3252  MOVE  L0G(.7),L0G(.2) 

3253  LABEL  "FROM" 

3254  MOVE  L0G(.7),L0G(.l) 

3255  LABEL  "EXPERIMENT" 

3500  LOIR  0 


APPENDIX  B 

OUTPUT  FROM  COMPUTER  PROGRAM  "CEDAR" 
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OUTPUT  FROM  COMPUTER  PROGRAM  "CEDAR" 


********************************************************************** 
★**★★*****★****★★*★★★*★*★**★★  CEDAR  *********>  **************** 
****************  CYLINDER  EXPANSION  DATA  REDUCTION  PROGRAM  ********* 
********************************************************************** 


*********  IMPUT  CYLINDER  EXPANSION  DATA  ********* 


I 

R-RO 

T 

Va 

Aa 

1 

0.0000 

-2.0000 

0.0000 

0.0000 

2 

BWmTTI 

3 

0.0000 

-.1400 

.7396 

.6688 

4 

1.0000 

1.2120 

1.0482 

.2676 

5 

2.0000 

2.1660 

1.2034 

.1738 

6 

3.0000 

2.9970 

1.2963 

.1159 

7 

4.0000 

3.7685 

1.3586 

.0827 

8 

5.0000 

4.5045 

1.4055 

.0648 

9 

6.0000 

5.2160 

1.4418 

.0517 

10 

7.0000 

5.9096 

1.4707 

.0421 

11 

8.0000 

6.5895 

1.4942 

.0348 

12 

9.0000 

7.2588 

1.5136 

.0292 

13 

10.0000 

7.9194 

1.5308 

.0262 

14 

11.0000 

8.5727 

1.5465 

.0240 

15 

12.0000 

9.2193 

1.5608 

.0223 

16 

13.0000 

9.8600 

1.5741 

.0208 

17 

14.0000 

10.4953 

1.5863 

.0193 

18 

15.0000 

11.1257 

1.5976 

.0180 

19  16.0000 

11.7516 

1.6077 

.0162 

20 

17.0000 

12.3736 

1.6168 

.0147 

21 

18.0000 

12.9921 

1.6254 

.0138 

22 

19.0000 

13.6074 

1.6334 

.0131 

23 

20.0000 

14.2196 

1.6410 

.0126 

24 

21.0000 

14.8290 

1.6484 

.0121 

25 

22.0000 

15.4356 

1.6555 

.0117 

26 

23.0000 

16.0397 

2.3882 

-.0949 

I  =■  INDEX  OF  DATA  POINT  READING  FROM  FILM  RECORD 
R-RO  »  MEASURED  EXTERIOR  WALL  DISLACEMENT  (MM) 

T  «  MEASURED  EXPANSION  TIME  (MICROSECONDS) 

Va  =  APPARENT  EXTERIOR  WALL  VELOCITY  (MM/MICROSEC) 

Aa  «  APPARENT  EXTERIOR  WALL  ACCELERATION  (W/MICROSEC®2) 


****************************  CALCULATED  DATA  **************************** 


I 

Rm-RmO 

Z 

Tp 

Vr 

Vz 

Ar 

Ac  Rho 

P 

Gaima 

2 

0.000 

0.000 

1.000 

0.000 

0.000 

0.000 

0.000  .00235 

29.50 

2.71 

3 

0.000 

0.000 

1.860 

0.000 

0.000 

0.000 

0.000  .00235 

29.50 

2.71 

4 

1.082 

.054 

3.219 

.844 

.042 

.621 

.031  .00184 

15.17 

2.72 

5 

2.154 

.129 

4.182 

1.128 

.080 

.295 

.039  .00141 

6.53 

3.13 

6 

3.217 

.214 

5.024 

1.271 

.102 

.170 

.026  .00116 

3.44 

3.32 

7 

4.273 

.304 

5.807 

1.354 

.116 

.105 

.018  .00099 

1.98 

3.43 

8 

5.323 

.397 

6.554 

1.407 

.125 

.072 

.013  .00086 

1.25 

3.26 

9 

6.368 

.492 

7.278 

1.447 

.132 

.054 

.010  .00075 

.89 

2.70 

10 

7.409 

.589 

7.983 

1.476 

.138 

.042 

.008  .00067 

.64 

2.75 

11 

8.446 

.688 

8.676 

1.499 

.142 

.033 

.006  .00060 

.48 

2.71 

12 

9.480 

.787 

9.357 

1.518 

.146 

.027 

.005  .00054 

.36 

2.70 

13 

10.511 

.887 

10.031 

1.532 

.149 

.022 

.004  .00049 

.28 

2.64 

14 

11.540 

.987 

10.696 

1.546 

.151 

.020 

.004  .00045 

.24 

1.79 

15 

12.566 

1.088 

11.356 

1.557 

.153 

.018 

.004  .00041 

.21 

1.52 

16 

13.590 

1.190 

12.009 

1.568 

.156 

.017 

.003  .00038 

.19 

1.42 

17 

14.613 

1.292 

12.657 

1.579 

.158 

.016 

.003  .00035 

.17 

1.41 

18 

15.634 

1.395 

13.300 

1.588 

.160 

.015 

.003  .00033 

.15 

1.55 

19 

16.654 

1.498 

13.939 

1.597 

.161 

.014 

.003  .00030 

.13 

1.45 

20 

17.673 

1.601 

14.574 

1.605 

.163 

.012 

.002  .00028 

.11 

2.35 

21 

18.690 

1.705 

15.206 

1.611 

.164 

.011 

.002  .00026 

.10 

2.15 

22 

19.706 

1.809 

15.834 

1.618 

.166 

.010 

.002  .00025 

.09 

1.34 

23 

20.722 

1.913 

16.459 

1.624 

.167 

.010 

.002  .00023 

.08 

1.45 

24 

21.736 

2.018 

17.082 

1.630 

.168 

.010 

.002  .00022 

.08 

1.06 

25 

22.750 

2.123 

17.702 

1.636 

.169 

.009 

.002  .00021 

.07 

1.03 

26 

23.763 

2.228 

18.319 

1.642 

.171 

.009 

.002  .00020 

.07 

.95 

Rm-RmO  »  RADIAL  DISPLACEMENT  OF  WALL  PARTICLE  CENTER  OF  MASS  (W) 
Z  »  AXIAL  DISPLACEMENT  OF  WALL  PARTICLE  CENTER  OF  MASS  (W) 
Tp  =  ELAPSED  TIME  SINCE  START  OF  PARTICLE  MOTION  (MICROSEC) 

Vr  =  LAGRAN6IAN  RADIAL  VELOCITY  COMPONENT  (MM/MICROSEC) 

Vz  =  LAGRANGIAN  AXIAL  VELOCITY  COMPONENT  (W/MICROSEC) 

Ar  *  LAGRANGIAN  RADIAL  ACCELERATION  COMPONENT  (MM/MICROSEC^) 
Az  *  LAGRANGIAN  AXIAL  ACCELERATION  COMPONENT  (MM/MICROSEC©2) 
Rho  *  DETONATION  PRODUCTS  GAS  DENSITY  (GRAMS/MM«3) 

P  =  DETONATION  PRODUCTS  GAS  PRESSURE  (GIGAPASCALS) 

Gairnia  =  GENERALIZED  RATIO  OF  SPECIFIC  HEATS 

END  OF  COMPUTATION 


Note:  The  symbol  "®"  denotes  exponentiation 


LIST  OF  SYMBOLS 


Symbol  Meaning  Page 

Aq  Cross-sectional  area  of  cylinder  wall  13 

D  Detonaton  velocity  of  explosive  13 

f(t)  Axial  displacement  17 

g(t)  Radial  displacement  17 

h(z,t)  Central  surface  displacement  observed  at  station  z,time  t  13 

£(t)  Arclength  along  path  of  particle  motion  17 

M  Mass  of  cylinder  per  unit  length  19 

P  Pressure  of  detonation  product  gases  20 

r  Eulerian  radial  coordinate  13 
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v  Radial  velocity  component  of  cylinder  wall  motion  22 
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z  Eulerian  axial  coordinate  13 

T  Generalized  ratio  of  specific  heats  21 

p  Density  of  detonation  product  gases  20 
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